/*******************************************************************************
                                HKLUTIL.C

  Title:        Hkl utility package for Speden.
  Author:       Hanna Szoke
  Date:         2/21/95 (original)
  Function:     This package contains functions that handle (hkl) factors 
		for all codes that need them except functions that read 
		and write structure factors.

 They include:
		low-level functions handling data -
			getAmpPhase()
			cleanPhase()
			fetch_hkl()
			fetch_Lhkl()
			assemble()
			fsetDsqLimit()
			scalefobs()
		R factor functions -
			get1Rfac()
			Lb_get1Rfac()
		Table Header functions -
			hklTabHead()
			shellHeader()
			shellHeader_R()
			shellHeader_F()
		Mask-related functions -
			get_P1mask()
			makeit_unique()
			prepare_unique_hklmask()
			check_others()		(internal to hklutil.c)
			check_friedel()		(internal to hklutil.c)
			set_centrics()
			set_cen_phases()
		and set-up functions formerly in hklinit.c-
			checkSigmas()
			SigmaRange()
			getNormFac()
			setupexpfac() 
			describeFR()
*******************************************************************************/

#include "speden.h"     /* includes util.h, cellparams.h, symmetry.h, dims.h */
#include	"sp.h"                 /* contains fobs, etc.    */

int	Nfob ;		/* length of Fo list	*/
int	Nfc ;		/* length of Fc list	*/
float	*Hlist ;	/* fractional H indices	for Fobs*/
float	*Klist ;	/* fractional K indices	for Fobs*/
float	*Llist ;	/* fractional L indices	for Fobs*/
real	*Folist ;	/* ungridded Fo values	*/
real	*Siglist ;	/* ungridded Sig values	*/
float	*Fc_Hlist ;	/* fractional H indices	for Fcalcs*/
float	*Fc_Klist ;	/* fractional K indices	for Fcalcs*/
float	*Fc_Llist ;	/* fractional L indices	for Fcalcs*/
real    *Fc_relist ;    /* ungridded Fc real values     */
real    *Fc_imlist ;    /* ungridded Fc imag. values    */

  /************************************************************
  declarations of EXTERN variables in speden.h               
  ************************************************************/

real  *hkl_weight;        	/* weights for (hkl) in half-ellipsoid */
real  *expfac ;       	/* exponential factors for sc solutions */

real	dsq_limit ;
static	char	Arange[NUMSEGS][MAXSTRING] ;
static	char	*unique_hkl ;
static  int	gridded_Nhkl ;

static	void	check_others(char *, int, int, int, int) ;
static	void	check_friedel(char *, int, int, int, int) ;

	/**************************************************
	Convert a complex (re, im) R into an amplitude and 
	phase; disallow very small amplitudes and put 
	phases into the expected range. 
	**************************************************/

void getAmpPhase(COMPLEX *R, real amp0, real *amplitude, real *phase) 
{
     real	cleanPhase(real *) ;
     real 	omega ;

     *amplitude = Cabs(R) ; 

     if (*amplitude < amp0*AMP_EPSI)  {
          *amplitude = 0 ;
          *phase = 0 ;
     }
     else { 
	  omega = RTOD * atan2(R->im, R->re) ;
          *phase = cleanPhase(&omega) ;
     }

}
real	cleanPhase(real *phase)
{
     while (*phase > 180.)
          *phase -= 360 ;
     while (*phase<= -180.)
          *phase += 360 ;
     if (fabs(180 - fabs(*phase)) < EPSI_PHASE)
          *phase = 180 ;
     else if (fabs(*phase) < EPSI_PHASE)
          *phase = 0 ;

     return(*phase) ;
}

	/**************************************************
	Rearrange negative k and l indices.
	Use position: 0, 1, 2, ... N-1,  N,  N+1, ... 2*N-1  
	for k (or l): 0, 1, 2, ... N-1, -N, -N+1, ... -1                
	**************************************************/

void fetch_hkl(int n, int *h, int *k, int *l) 
{

          *h = n % Nh ;
     	  *k = (n/Nh) % Nk ;
     	  *l = n / (Nh*Nk) ;

	  if (*k >= Nk/2)
	       *k -= Nk ;

	  if (*l >= Nl/2)
	       *l -= Nl ;

}
void fetch_Lhkl(int n, float *h, float *k, float *l)
{
     *h = *(Fc_Hlist + n) ;
     *k = *(Fc_Klist + n) ;
     *l = *(Fc_Llist + n) ;
}
int	assemble(int h, int k, int l) 
{

      if (h < 0 || h >= Nh)
          return (-1) ;
      if ((k < -Nk/2) || (k >= Nk/2))
          return (-1) ;
      if ((l < -Nl/2) || (l >= Nl/2))
          return (-1) ;

      return (h + ((Nk + k) % Nk)*Nh + ((Nl + l) % Nl)*Nh*Nk) ;
}

void	fsetDsqLimit(int N, char fo)
{
     real	fo_nusq(int) ;
     real	fc_nusq(int) ;
     real	nusqval = 0 ;
     static	char	DONE_IT = FALSE ;
     int	n ;
     real	thisnusq ;

     if (!DONE_IT) {

          /* Fetch max(nusq) for N entries.  */

        if (fo) 
          for (n = 0; n < N; n++) {
               thisnusq = fo_nusq(n) ;
	       if (thisnusq > nusqval) 
	           nusqval = thisnusq ; 
          }
        else    
          for (n = 0; n < N; n++) {
               thisnusq = fc_nusq(n) ;
	       if (thisnusq > nusqval) 
	           nusqval = thisnusq ; 
          }
          fprintf(fp_log, "\nNusqval = %g\n", nusqval) ;
          if (nusqval < limit_rat) {
	       dsq_limit = nusqval ;
               fprintf(fp_log, 
	  "\n1/d-squared shell limit is %g, based on data\n\n", eta*dsq_limit) ;
          }
          else {
	       dsq_limit = limit_rat ;
               fprintf(fp_log, 
	  "\n1/d-squared shell limit is %g, based on Speden parameter\n\n", 
	       eta*dsq_limit) ;
          }
          DONE_IT = TRUE ;
     }
}

		/*********************************************
	        Calculate the R factor segment by segment.
		R factor = Sum [ |Fob - |R|| ] / Sum [ Fob ]
		*********************************************/

real get1Rfac(real *Flocal, real *Rlocal, int Nlocal)
{
     real	fo_nusq(int) ;
     real	sumFR[NUMSEGS] ;
     real	sumF[NUMSEGS]  ;
     real	Rabs ;
     real	Fabs ;
     real	rat[NUMSEGS] ;
     real	totFR = 0, totF = 0, totrat ;

     int	m, n ;

     for (m = 0; m < NUMSEGS; m++) {
	  sumFR[m] = 0 ;
	  sumF[m] = 0 ;
     }

     for (n = 0; n < Nlocal; n++) {

          for (m = 1; m < NUMSEGS; m++)
	    if (fo_nusq(n) < nusqlim[m])
		 break ;

	  if (m < NUMSEGS) {   

               Rabs = fabs(*(Rlocal+n)) ;
	       Fabs = fabs(*(Flocal+n)) ;

	       sumFR[m] += fabs(Fabs - Rabs) ;
               sumF[m]  += Fabs ;
	  }
     }
     /*** NOTE: Ther is NO (000) value for List input!  ***/

     for (m = 1; m < NUMSEGS; m++) {
	  totFR += sumFR[m];
	  totF  += sumF[m];
	  if (sumF[m] != 0)
	       rat[m] = sumFR[m] / sumF[m] ;
          else
	       rat[m] = 0. ;
     }

     totrat = totFR / totF ;

     sprintf(message, "%s %10g %10g %10g %10g %10g", 
	    "      ", rat[1], rat[2], rat[3], rat[4], totrat) ;
     printTwice(message) ;

     return (totrat) ;
}

void	scalefobs(real *list, real sfac)
{
     int	m ;

     if (sfac != 1.) 
	  for (m = 0; m < Nfob; m++, list++) 
	       *list *= sfac ;

}
void hklTabHead()
{
     real	Amin ;
     int	m ;

     Amin = sqrt(1 / dsq_limit) ;

	 /* Initialize segment limits in Angstrom with "nice" numbers. */

     for (m = 1; m < NUMSEGS; m++) 
          sprintf(Arange[m], "%4.1f", pow(2.0, 0.5*(NUMSEGS-m-1))*Amin) ;

   /* Recalculate limits in terms of 1/d^2, but don't fiddle with outermost! */

     for (m = 1; m < NUMSEGS; m++) {	 

          sscanf(Arange[m], "%g", &t1) ;
          nusqlim[m] = 1. / (t1*t1) ;

	  if (m == NUMSEGS-1)		
               nusqlim[m] = dsq_limit ;	

     }
}
		/*********************************************
	        Calculate the R factor segment by segment.
		using lists of off-grid F's and R's. 
		*********************************************/

real Lb_get1Rfac(COMPLEX *Oar)
{
     real	fc_nusq(int) ;
     COMPLEX	cdiff ;
     real	modsq, diffsq ;
     real	resq, imsq ;
     real	sumFR[NUMSEGS] ;
     real	sumF[NUMSEGS]  ;
     real	rat[NUMSEGS] ;
     real	totFR = 0, totF = 0, totrat ;
int totc = 0 ;
int cases[NUMSEGS] ;
     int	m, n ;

     for (m = 0; m < NUMSEGS; m++) {
	  sumFR[m] = 0 ;
	  sumF[m] = 0 ;
cases[m] = 0 ;
     }

     for (n = 0; n < Nfc; n++) {

	  if (n == 0)
               m = 0 ;

          else {
	       for (m = 1; m < NUMSEGS; m++)
		    if (fc_nusq(n) < nusqlim[m])
			 break ;
          }
	  if (m < NUMSEGS) {   

                cdiff.re =  (*(Fc_relist+n) - (Oar+n)->re) ;
                cdiff.im =  (*(Fc_imlist+n) - (Oar+n)->im) ;
                diffsq = cdiff.re * cdiff.re + cdiff.im * cdiff.im ;

                resq = (*(Fc_relist+n) * *(Fc_relist+n)) ;
                imsq = (*(Fc_imlist+n) * *(Fc_imlist+n)) ;
	        modsq = resq + imsq ;

	        sumFR[m] += sqrt(diffsq) ; 
                sumF[m]  += sqrt(modsq) ;
cases[m]++ ;
	  }
     }

     for (m = 0; m < NUMSEGS; m++) {
	  totFR += sumFR[m];
	  totF  += sumF[m];
totc += cases[m] ;
          if (sumF[m] > 0)
	       rat[m] = sumFR[m] / sumF[m] ;
          else
	       rat[m] = 0 ;
     }

     if( totF != 0 )
	totrat = totFR / totF ;
     else 
	totrat = 0.;

     sprintf(message, "Model:  %10g %10g %10g %10g %10g %10g", 
	    rat[0], rat[1], rat[2], rat[3], rat[4], totrat) ;
     printTwice(message) ;

     sprintf(message, "# cases:%10d %10d %10d %10d %10d %10d",
	    cases[0], cases[1], cases[2], cases[3], cases[4], totc) ;
     printTwice(message) ;

     return (totrat) ;
}

void shellHeader()	/* called by dphase and solve; no (000) */
{
     printTwice("\n\t\tR factors for fractions of 1/d-squared\n") ;
     printTwice(
     "\nhkl shell   <1/8    1/8-1/4      1/4-1/2       rest") ;

     sprintf(message,
     "resol (nm)  %s  %s-%s  %s-%s  %s-%s    all \n",
     Arange[1], Arange[1], Arange[2], Arange[2], Arange[3], Arange[3],
     Arange[4]) ;
     printTwice(message) ;
}
void shellHeader_R()	/* called by Back;  */
{
     printTwice("\n\t\tR factors for fractions of 1/d-squared\n") ;
     printTwice(
    "\t       (000)        1/8        1/4        1/2  Remainder    Overall\n") ;

     sprintf(message, 
     "resol (nm)               > %s  %s-%s  %s-%s  %s-%s        all \n",
     Arange[1], Arange[1], Arange[2], Arange[2], Arange[3], Arange[3],
     Arange[4]) ;
     printTwice(message) ;
}

void shellHeader_F()	/* called by init (for Solve only); goes to log. */
{
     fprintf(fp_log,
     "           (000)   1/8        1/4        1/2  Remainder    Overall\n") ;
     fprintf(fp_log, 
     "res (nm)          > %s  %s-%s  %s-%s  %s-%s        all \n\n",
     Arange[1], Arange[1], Arange[2], Arange[2], Arange[3], Arange[3],
     Arange[4]) ;
}
	/*********************************************
        Create a mask in which we enforce the
	half-ellipsoid rule (for expanded output).
	*********************************************/

char	*get_P1mask() 
{
     char	*mask ;
     int	n, h, k, l ; 
     gridded_Nhkl = Nh*Nk*Nl ;

     mask = (char *) e_malloc(gridded_Nhkl*sizeof(char), "get_P1mask") ; 
          
     for (n = 0; n < gridded_Nhkl; n++) {

          *(mask+n) = 0 ;

          fetch_hkl(n, &h, &k, &l) ;

          if (nusq(h, k, l) < limit_rat) {
	  
	       if (((h > 0)) ||
		   ((h == 0) && (k > 0)) ||
		   ((h == 0) && (k == 0) && (l >=0))) 
	            *(mask+n) = 1 ;
          }
     }

     return (mask) ;
}

	/*********************************************
        Check a mask, resetting values that are not 
	part of the unique set.  Return the number of
	reset values, if changes were made.
	*********************************************/

int	makeit_unique(char *mask) 
{
     void	prepare_unique_hklmask() ;
     int	n ;
     char	*unique_mask ;
     int	nonu = 0;
     int	u = 0 ;

     prepare_unique_hklmask(&unique_mask) ;

     for (n = 1; n < Nfc; n++)
	  u += *(mask+n) ;

     for (n = 1; n < Nfc; n++)

	  if (*(mask+n) && !*(unique_mask+n)){
	       nonu++ ;
	       *(mask+n) = 0 ;
          }

     e_free(Nfc*sizeof(char), unique_mask) ;

     return (nonu) ;

}

		/**************************************************
		Determine unique reflections, for which **pmask = 1.  
		Procedure: Initialize all points in mask to zero.
		Set unique points to 1, related non-unique pts to 2.
		When all mask has been accounted for, reset to 0 
		all points outside the resolution limits and
		all points that were previously set to 2.
		**************************************************/

void prepare_unique_hklmask(char **pmask)
{
     int	n ;
     float	h, k, l ;
     int	kk, ll ;
     int	sum_as1 = 0, sum_notas1 = 0 ;
     char	*umask ;

     unique_hkl = (char *) e_malloc(Nfc*sizeof(char), "hklutil") ;
     
     *pmask = unique_hkl ;
     
     for (n = 0, umask = unique_hkl; n < Nfc; n++, umask++)
	  *umask = 0 ;

     umask = unique_hkl ;

     for (h = 0; h < Nh; h++) {

          for (kk = 0; kk < Nk; kk++) {

	       k =  (kk > Nk/2) ? kk - Nk : kk ;

	       if (k != Nk/2) {

                    for (ll = 0; ll < Nl; ll++) {

	                 l =  (ll > Nl/2) ? ll - Nl : ll ;

		         if (l != Nl/2) {

                              if ((n = assemble(h, k, l)) >= 0) {

			           if (*(umask+n) == 1) {
                              sprintf(message, 
	   		      "Illegally set point at hkl = %f %f %f",
     			      h, k, l) ;
                              Fatal(message) ;
			           }
			           else if (*(umask+n) == 0) {
			                *(umask+n) = 1 ;
			                check_others(unique_hkl, h, k, l, n) ;
                                   }
			           /*  else if (*(umask+n) == 2) do nothing */
                              }
                         }
                    }		
               }		
          }	
     }	
     
     /* Reset the mask outside the nusq limit */

     for (n = 0, umask = unique_hkl; n < Nfc; n++, umask++) {

          if (fc_nusq(n) >= limit_rat)  
	       *umask = 0 ;

     }			

     /* Statistics */

     for (n = 0, umask = unique_hkl; n < Nfc; n++, umask++) {
	  if (*umask == 1)
	       sum_as1++ ;
	  else if (*umask == 2)
	       sum_notas1++ ;
	  /* else (*umask == 0) not interested */
     }

     sprintf(message, "\nMax. # of unique and non-unique reflections: %d %d\n", 
	     sum_as1, sum_notas1) ;
     printTwice(message) ;


     /* Now replace "others" by "masked-out" */

     for (n = 0, umask = unique_hkl; n < Nfc; n++, umask++) 
	  if (*umask == 2)
	       *umask = 0 ;
}

void	check_others(char *umask, int h, int k, int l, int n)
{
     int	p ;
     real	off ;
     real	*mat ;
     int	newh, newk, newl, newn ;

     for (p = 0, mat = matop; p < Ncs; p++, mat += MEL) {

          apply_symop_hkl(h, k, l, mat, &newh, &newk, &newl, &off) ;

          if (((newn = assemble( newh,  newk,  newl)) >= 0) &&
               (newn != n)) {

		if (*(umask+newn) == 1) {
                     sprintf(message, 
     "Illegally set unique point - hkl = %d %d %d, new h,k,l = %d %d %d", 
                     h, k, l, newh, newk, newl) ;
                     Fatal(message) ;
                }
		else if (*(umask+newn) == 0) 
		     *(umask+newn) = 2 ;
		/* else if (*(umask+newn) == 2) do nothing */
          }
          check_friedel(umask, newh, newk, newl, n) ;
     }
}

void	check_friedel(char *umask, int h, int k, int l, int oldn)
{
     int	n ;

     if (((n = assemble(-h, -k, -l)) > 0) && (n != oldn)) {

         if (*(umask+n) == 1) {
              sprintf(message, "Illegally set Friedel point at hkl = %d %d %d",
		           -h, -k, -l) ;
               Fatal(message) ;
         }
         else if (*(umask+n) == 0) {
              *(umask+n) = 2 ;
         }
         /*  else if (*(umask+n) == 2) do nothing */
     }
}

void	set_centrics(char *mcentric)
{
     real	off ;
     real	*mat ;
     int	N ;
     int        p, n ;
     float      h, k, l ;
     int	newh, newk, newl ;

     for (n = 0; n < Nfc; n++) 
          *(mcentric+n) = 0 ;

	/***********************************************************
	We look for (hkl)'s which, under a symmetry operation, 
	give back the negatives of the same reciprocal indices.
	See Giacovazzo, Fundamentals of Crystallography, pp 156-157. 
	***********************************************************/

     for (n = 1; n < Nfc; n++) {

          fetch_Lhkl(n, &h, &k, &l) ;

          if (fc_nusq(n) < limit_rat) {

               for (p = 0, mat = matop; p < Ncs; p++, mat += MEL) {

	           apply_symop_hkl(h, k, l, mat, &newh, &newk, &newl, &off) ;

	           if ((newh == -h) && (newk == -k) && (newl == -l)) {
		       *(mcentric+n) = 1 ;
		       break ;
                   }		
               }	
          }	
     }			

     for (n = 0, N = 0; n < Nfc; n++) 
          N += *(mcentric+n) ;

     /*** return(N) ; ***/
}

int set_cen_phases(char *mcentric, real *mphases)
{
     real	off ;
     real	*mat ;
     int	N ;
     int        p, n ;
     float      h, k, l ;
     int	newh, newk, newl ;

     for (n = 0; n < Nfc; n++) {
          *(mcentric+n) = 0 ;
          *(mphases+n) = 0 ;
     }

	/***********************************************************
	We look for (hkl)'s which, under a symmetry operation, 
	give back the negatives of the same reciprocal indices.
	We set the corresponding phases according to Giacovazzo, 
	3.39 (but without the n*pi factor that is added in ranphase)
	***********************************************************/

     for (n = 1; n < Nfc; n++) {

          fetch_Lhkl(n, &h, &k, &l) ;

          if (fc_nusq(n) < limit_rat) {

               for (p = 0, mat = matop; p < Ncs; p++, mat += MEL) {

	           apply_symop_hkl(h, k, l, mat, &newh, &newk, &newl, &off) ;

	           if ((newh == -h) && (newk == -k) && (newl == -l)) {
		       *(mcentric+n) = 1 ;
		       *(mphases+n) = off / 2. ;
		       break ;
                   }		
               }	
          }	
     }			

     for (n = 0, N = 0; n < Nfc; n++) 
          N += *(mcentric+n) ;

     return(N) ;
}

  	/************************************************************
	Examine sigmas:
	Reset useSig if all sigmas are 0. 
	Exit with error if some sigmas are set and others aren't.
	Then, if useSig is off, replace sigma = 0 by 1  everywhere.
	************************************************************/

int checkSigmas(real *sarray, int N)
{
     int	s ;
     real	*sig ;
     int	Nsigpresent = 0, Nsigabsent = 0 ;

     if (useSig) {

          for (s = 0, sig = sarray; s < N; s++, sig++)  {
	       if (*sig == 0) 
	            Nsigabsent++ ;   
	       else
		    Nsigpresent++ ;
          }
          if (Nsigpresent == 0) {
                printTwice("\nAll sigmas are 0 - turning off useSig flag") ;
		useSig = FALSE ;
          }
	  else if (Nsigabsent > 0) {
	       Fatal(
		     "Bad fobs file: some sigmas are 0, others are not.") ;
          }
     }

     if (!useSig) 
          for (s = 0, sig = sarray; s < N; s++, sig++)  
	       *sig = 1. ;
     
     return(useSig) ;

}

  	/************************************************************
	Determine range of sigmas; 
	Ignore sigma(000) for purposes of reporting range,
	5/1/01: there was an error and we had NOT ignored sigma(000)
	but now that error is corrected.
	************************************************************/

void	SigmaRange(real *sigmin, real *sigmax, real *sarray, int N)
{
     int	s ;
     real	smin = 1000. ;
     real	smax = -1000. ;

     for (s = 0; s < N; s++, sarray++) {
          if (*sarray < smin)
	      smin = *sarray ;
	  if (*sarray > smax)
	      smax = *sarray ;
     }          

     if (smin > 0)  {	/* legal range of sigmas */
          sprintf(message, "Sigma range is %g - %g\n", smin, smax) ;
          printTwice(message) ;
     }
     else {				/* something's wrong */
          sprintf(message, 
           "The range of sigmas: (%g, %g) is illegal!", smin, smax) ;
          Fatal(message) ;
     }

     *sigmin = smin ;
     *sigmax = smax ;

}

  	/************************************************************
        A single data set is used for normalization.
	As of 9/28/98. sigmas = 1 if useSig is not set.
	************************************************************/

real getNormFac() 	
{
     real	SsqWts ;	
     real	Smasks ;	
     int	n ;
     real	*wt ;
     real	*sig ;

     wt = hkl_weight ;
     sig = Siglist ;

     for (n = 0, SsqWts = 0, Smasks = 0; n < Nfob; n++, wt++, sig++) {
          SsqWts += (*wt / *sig) * (*wt / *sig) ; 
          Smasks += 1 ;
     }

     return (SsqWts / Smasks) ;
}

  	/************************************************************
        Set up expfac(N), defined as: exp (-A|k*dr|^2) where
        A is PISQ*eta, k stands for (h,k,l), dr for (dx, dy, dz). 
	dfac is	usually = delfac = PISQ * eta * dr * dr
	************************************************************/

real	*setupexpfac(real dfac)
{
     real	fc_nusq(int) ;
     real	*ef ;
     int        n ;

     strcpy(message, "setupexpfac") ;

     ef = (real *) e_malloc(Nfc*sizeof(real), message) ;  

     for (n = 0; n < Nfc; n++) 
          *(ef+n) = exp(-(dfac * fc_nusq(n))) ;

     return (ef) ;
}

real	*setupexpfac_fo(real dfac)
{
     real	fo_nusq(int) ;
     real	*efo ;
     int        n ;

     strcpy(message, "setupexpfac_fo") ;

     efo = (real *) e_malloc(Nfob*sizeof(real), message) ;  

     for (n = 0; n < Nfob; n++) 
          *(efo+n) = exp(-(dfac * fo_nusq(n))) ;

     return (efo) ;
}

	/*********************************************/
       	/* Calculate F and counts segment by segment.*/	
	/*********************************************/

void 	describeFR(char *maskc)
{
     real       fo_nusq(int n) ;
     void	shellHeader_F() ;
     void	hklTabHead() ;

     real	sumF[NUMSEGS+1] ;
     int	ocases[NUMSEGS+1], ccases[NUMSEGS+1] ;
     int	valid_cases[NUMSEGS+1] ;
     int	m, n ;
     int	h, k, l ;

     for (m = 0; m <= NUMSEGS; m++) {
	  sumF[m] = 0 ;
	  ocases[m] = 0 ;
	  ccases[m] = 0 ;
	  valid_cases[m] = 0 ;
     }
     hklTabHead() ;

     for (n = 0; n < Nfob; n++) {

	  for (m = 1; m < NUMSEGS; m++)
	     if (fo_nusq(n) < nusqlim[m])
	         break ;
          
	  if (m < NUMSEGS) {   
               sumF[m]  += *(Folist+n) ;
	       ocases[m] += 1 ;
	       valid_cases[m] += 1 ;
          }
     }

     for (n = 0; n < gridded_Nhkl; n++) {

          fetch_hkl(n, &h, &k, &l) ;

	  for (m = 1; m < NUMSEGS; m++)
	       if (nusq(h, k, l) < nusqlim[m])
		   break ;
          
	  if (m < NUMSEGS) {   
	       ccases[m] += *(maskc+n) ;
          }
     }

     sumF[0] = F000 ;
     ocases[0] = ccases[0] = valid_cases[0] = 1 ;

     for (m = 0; m < NUMSEGS; m++) {
	  sumF[NUMSEGS] += sumF[m];
	  ocases[NUMSEGS] += ocases[m];
	  ccases[NUMSEGS] += ccases[m];
	  valid_cases[NUMSEGS] += valid_cases[m];
     }

     for (m = 0; m <= NUMSEGS; m++) {
	  sumF[m] = sumF[m] * 100 / sumF[NUMSEGS] ;
     }

     fprintf(fp_log,
     "\nThe following table should help you estimate to what extent\n") ;
     fprintf(fp_log,
     "your problem is under- or overdetermined.  For each resolution bin,\n") ;
     fprintf(fp_log,
     "numbers in row # Fo count the entries in your input Fobs;\n") ;
     fprintf(fp_log,
     "numbers in row # Fc count the integer grid points that would be\n") ;
     fprintf(fp_log,
     "in a fully-gridded fcalc file at the given resolution.\n\n") ;

     shellHeader_F() ;

     fprintf(fp_log, 
            "%% F %9g %10g %10g %10g %10g %10g\n", 
	    sumF[0], sumF[1], sumF[2], sumF[3], sumF[4], sumF[5]) ;
     fprintf(fp_log, 
            "# Fo    %5d %10d %10d %10d %10d %10d\n", 
	    ocases[0], ocases[1], ocases[2], ocases[3], ocases[4], ocases[5]) ;
     fprintf(fp_log, 
            "# Fc    %5d %10d %10d %10d %10d %10d\n", 
	    ccases[0], ccases[1], ccases[2], ccases[3], ccases[4], ccases[5]) ;

     fflush(fp_log) ;
     return ;
}
